#pragma once
#include "Vector.hpp"
#include "Matrix.hpp"
#include "../common.hpp"
#include "Math.hpp"
#include "SymEigenSystem.hpp"
#include "SumCount.hpp"

namespace zzz{
template<typename T>
T Mean(const vector<T> &data)
{
  SumCount<T> sc;
  sc += data;
  return sc.Mean();
}

template<typename T>
T Variance(const vector<T> &data)
{
  T mean=Mean(data);
  T var(0);
  for (zuint i=0; i<data.size(); i++) {
    T diff=data[i]-mean;
    var+=diff*diff;
  }
  return var/(data.size()-1);
}

template<typename T>
T StdDeviation(const vector<T> &data)
{
  return Sqrt<T>(Variance(data));
}

template<zuint N, typename T>
Vector<N, T> Covariance(const vector<Vector<N, T> > &data, MatrixBase<N, N, T> &cov)
{
  Vector<N, T> mean=Mean(data);
  vector<Vector<N, T> > diff;
  diff.reserve(data.size());
  for (zuint i=0; i<data.size(); i++)
    diff.push_back(data[i]-mean);
  cov.Zero();
  for (zuint i=0; i<diff.size(); i++) {
    for (zuint n=0; n<N; n++) for (zuint m=n; m<N; m++)
      cov(n, m)+=diff[i][n]*diff[i][m];
  }
  cov/=(data.size()-1);
  for (zuint n=1; n<N; n++) for (zuint m=0; m<n; m++)
    cov(n, m)=cov(m, n);
  return mean;
}

template<zuint N, typename T>
Vector<N, T> PCA(const vector<Vector<N, T> > &data, MatrixBase<N, N, T> &eigen_vector, Vector<N, T> &eigen_value)
{
  Matrix<N, N, T> cov;
  Vector<N, T> mean = Covariance(data, cov);
  SymEigenSystem<T, N> ss;

  ss.Decompose(cov, eigen_vector, eigen_value);
  return mean;
}
}
